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1.  INTRODUCTION:  A  BIT  OF  HISTORY 

"...no  one  doubts  that  the  modern  formulations  (of  science)  are  clear, 

elegant,  and  precise;  it's  just  that  it's  impossible  to  comprehend  how 

anyone  ever  thought  of  them." 

— Michael  Spivak,  A  Comprehensive  Introduction  to  Differential  Geometry 

In  the  last  few  years,  the  word  "fractal"  (a  set  with  noninteger 
fractional  dimension)  has  worked  its  way  out  from  research  into  more  general 
use  in  the  scientific  community.  The  word  was  coined  in  the  seventies  by 
Benoit  B.  Mandelbrot  to  describe  the  elaborate  images  he  was  producing  on 
comr  ,;ers  at  the  Watson  IBM  Research  Center.  However,  as  Mandelbrot  himself 
poinv=>  out  in  his  book,  The  Fractal  Geometry  of  Nature  (1983),*  the  roots  of 
the  fractal  idea  reach  back  over  a  century. 

In  the  nineteenth  century,  mathematicians  were  working  to  develop  rigor  in 
their  study  of  mathematics.  They  encountered  curious  phenomena,  which  forced 
the  rethinking  of  some  concepts.  Four  of  these  phenomena  are  cited: 

1.  In  1874,  K.  Weierstrass  produced  a  continuous,  but  nowhere  differen¬ 
tiable  function.  To  create  this  function,  he  employed  trigonometric  series 
and  lacunary  (or  "gap")  series.  An  example  of  a  Weierstrass  function  (fig.  1) 
is  given  by  the  formula 

f(t)  =  l  [(|)n]  cos ( 2nt )  . 
n 

(Weierstrass'  original  results  were  more  general.) 

2.  In  the  l870's,  0.  Cantor  produced  several  results  which  gave  relation¬ 
ships  between  set  theory  and  calculus.  In  1 87 ^ ,  he  produced  his  proof  that 
there  are  only  countably  many  algebraic  numbers,  and  became  increasingly 
fascinated  with  the  concept  of  infinity  in  mathematics.  A  concrete  example  of 
his  ideas  is  the  "middle  thirds"  set,  which  was  used  by  H.  Lebesgue  in  the 
1920's  to  produce  the  Cantor-Lebesgue  function  (fig.  2).  This  singular  func¬ 
tion  has  a  derivative  equal  to  zero  at  almost  all  points  in  the  set  (0,1),  yet 
is  monotonically  increasing. 

3.  In  1890,  G.  Peano  produced  an  example  of  a  continuous  space-filling 
curve  (fig.  3).  This  curve  maps  the  unit  interval  [0,1]  onto  the  unit  square 
[0,1]  *  [0,1]. 

4.  In  1904,  h.  von  Koch  produced  his  snowflake  (fig.  4).  This  curve  has 
infinite  length,  but  is  contained  in  a  finite  area.  The  snowflake  is  non- 
intersecting  and  is  self-similar,  i.e.,  it  appears  the  same  despite  successive 
magnifications. 

In  each  of  the  above  examples,  a  rigorous  limit  procedure  was  used  (uni¬ 
form  convergence  in  the  proper  topology).  However,  because  each  of  these 
examples  conflicted  with  the  geometric  intuition  of  the  time,  they  were 


*See  bibliography. 


axis 


labeled  "pathological."  Yet,  because  of  further  scientific  research,  and 
especially  because  of  the  introduction  of  the  computer,  it  is  now  possible  to 
see  that  these  examples  appear  to  model  nature  quite  well.  Each  curve  fits 
Mandelbrot’s  basic  definition  of  a  fractal  curve:  a  curve  having  fractional 
dimension  higher  than  one.  Mandelbrot  cites  all  of  them  as  important  examples 
in  his  book. 

In  the  following  sections,  theory  and  procedure  for  the  creation  of  frac¬ 
tal  images  are  discussed.  The  discussion  is  aimed  at  the  reader  with  some 
computer  graphics  experience  and/or  software.  Numerous  fractal  images  are 
included,  and  algorithms  for  the  production  of  some  of  these  images  are  pro¬ 
vided.  The  algorithms  can  be  set  up  to  run  on  almost  any  system — whenever 
possible,  hardware  specifics  are  eliminated.  Thus,  readers  can  run  these 
algorithms  on  their  own  systems  and  proceed  to  explore  the  world  of  fractals. 
Discussion  is  divided  into  three  areas:  fractals  on  complex  dynamical  sys¬ 
tems,  seed  fractals,  and  squig  fractals. 

The  theoretical  discussions  in  the  following  sections  are  independent  of 
the  descriptions  of  algorithms  which  produce  fractal  images.  In  particular, 
theoretical  discussions  precede  algorithms  in  section  2  (on  complex  analytic 
dynamics).  These  were  provided  to  give  the  interested  reader  some  insight  as 
to  why  the  algorithms  work.  Also,  section  3.2  (on  dimension)  must  be  labeled 
"theoretical,"  as  discussion  and  calculation  here  are  rather  involved.  Theo¬ 
retical  sections  can  be  skipped  over  with  a  minimal  loss  of  continuity. 

2.  COMPLEX  ANALYTIC  DYNAMICS 

Many  pioneers  besides  those  already  mentioned  were  involved  in  the  evolu¬ 
tion  of  the  fractal--Riemann,  Hausdorff,  Klein,  Cesaro,  and  Bernoulli,  to  name 
a  few — the  list  is  very  long.  One  research  field  related  to  fractals  that  has 
had  quite  a  revival  lately  is  the  field  of  iteration  theory,  or  complex  analy¬ 
tic  dynamics.  This  was  started  in  the  early  1900's  by  P.  Fatou  and  G.  Julia. 
Both  wrote  long  monographs  on  the  subject.  Today,  its  pioneers  include  D. 
Sullivan,  J.  Hubbard,  and  A.  Douady. 

2.1  Theory 

A  dynamical  system  consists  of  a  pair,  (X,  4>) ,  where  X  is  a  topologi¬ 
cal  space,*  and  4>  =  { <f>t :  t  €  R}  is  a  set  of  dynamics,  i.e.,  rules  for  the 
evolution  of  the  system  in  time.  If  4>t  is  continuous,  the  pair  (X,  <j>t)  is 
usually  called  a  flow. 

Examples  of  dynamical  systems  appear  in  numerous  places:  the  cardio¬ 
vascular  system,  the  capitalist  system,  and  the  solar  system  are  all  dynamic. 
In  each,  a  process  is  occurring  which  can  be  thought  of  as  evolving  in  time. 


M  topological  space  is  a  set  in  which  the  concept  "open  set"  is  well  defined.  These  topics  are  discussed  by 
Munkres  { 1975),  bibliography. 


In  this^section,  the  underlying  space  of  the  dynamical  system  is  the 
Riemann  sphere,  C,  where 


C  =  C  U  {*}  =  (z=x+iy:  x ,  y  €  R ,  i  =  /-T}  U 


The  point  "*>"  is  added  to  C,  the  complex  plane,  by  rolling  the  plane  up  into  a 
sphere,  and  letting  *  be  the  north  pole.  The  dynamic  is  an  analytic  func¬ 
tion.*  Let  f(z)  denote  this  function.  Then,  by  iterating  the  function 


zn+1  '  f(zn)  * 


one  gets  a  dynamical  system.  Note  that  this  system  is  discrete. 


The  computer  has  proven  to  be  a  most  useful  tool  in  the  study  of  nonlinear 
dynamical  systems.  This  study  has  produced  as  a  byproduct  some  of  the  fractal 
images  seen  here.  Of  these,  many  have  been  produced  by  complex  analytic 
dynamics.  Probably  the  most  common  dynamic  in  generating  these  images  has 
been  the  now-famous  equation 


f(z)  =  z2  +  A 


Using  this  equation,  we  can  discuss  the  two  different  types  of  images  it 
produces.  The  first  type  is  a  C-dynamical  system.  In  this  process,  the 
number  A  is  kept  constant,  and  z  is  varied.  The  second  type  of  image  is  a 
parameter  space  image.  Here,  z  is  fixed,  and  the  number  A^  is  varied.  Each 
different  value  of  A  parameterizes  a  dynamical  system  from  C.  The  Mandelbrot 
set,  which  is  generated  in  this  fashion  using  z2  +  A,  is  a  parameter  space 
image. 


The  following  discussion  of  some  of  the  theory  behind  generating  these 
images  comes  under  the  category  "complex  quadratic  dynamics."  Blanchard 
(1984)  gives  a  more  complete  discussion  (see  bibliography,  "Mathematics,"  for 
other  necessary  background). 


First,  consider  the  dynamics  in  C.  Heuristically ,  the  Riemann  sphere, 
instead  of  just  the  complex  plane,  Ck  is  the  base  space  because  of  the  impor¬ 
tance  of  infinity  in  dynamics.  In  C,  the  group  of  one-on-one  analytic  map¬ 
pings  is  the  group  of  Moebius  transforms,  that  is,  maps  of  the  form 


g(z)  = 


az  +  b 


cz  +  d 


ad  -  be  *  0 


Moebius  transformat  ions  have  one  zero  and  one  pole.  Using  these  maps,  it  is 
possible  to  get  some  insight  into  why  the  single  equation  f(z)  =  z2  +  A  is  so 
powerful.  Let  h(z)  =  Az2  +  2Bz  +  C  be  a  general  quadratic  equation  in  C.  For 
f(z)  as  above,  it  is  possible  to  solve  for  a  Moebius  transform  g(z)  and  a 
value  of  A  in  f(z)  such  that 


h(z)  =  g_1  0  f  0  g(z)  =  g-1  (f (g(z)) ) 


M  function  is  called  analytic  in  some  open  set  if  it  can  be  expanded  in  a  Taylor  series  in  that  set.  In 
C-analysis,  if  a  function  has  a  continuous  derivative  in  an  open  set,  it  is  analytic  there. 


WWW3 


The  solution  is  given  by 


g(z)  =  Az  +  B  , 

X  =  AC  +  B  -  B2  . 

Since  two  dynamical  systems  in  C  conjugate  by  a  Moebius  transform  are  the 
same,  every  quadratic  dynamical  system  in  C  can  be  obtained  by  varying  X.  In 
other  words,  every  quadratic  system  is  parameterized  by  the  complex  number  X. 

Given  the  dynamic  f(z)  =  z2  +  X  for  fixed  X,  iteration  produces  a  dynam¬ 
ical  system.  Under  this  iteration,  individual  points  will  have  neighborhoods 
(small  open  sets  containing  the  point  in  discussion)  that  exhibit  one  of  two 
types  of  behavior.  Points  in  these  neighborhoods  will  either  converge  to  a 
point  after  repeated  iteration,  or  they  will  not  converge.  Those  points  that 
have  a  neighborhood  of  points  that  converge  are  called  elements  of  the  Fatou 
set.  The  points  not  in  the  Fatou  set  are  called  elements  of  the  Julia  set. 

The  function  f(z)  =  z2  provides  an  enlightening  example  of  this  dichot¬ 
omy.  Note  that  under  iteration  of  f,  every  point  with  absolute  value  strictly 
less  than  one  will  converge  to  zero,  while  every  point  z  with  |z|  >  1  will 
converge  to  infinity.  However,  under 

fn  =  f  °  ...  c  f  =  f (f ( . .  .  (f (z) ) ) )  ,  n  times, 

most  points  on  the  unit  circle  (i|z|  =  1})  are  just  "spun  around"  at  a  faster 
and  faster  rate.  In  fact,  if  z  =  exp(ia),  where  a  is  an  irrational  number, 
there  exists  an  iterate  of  z  coming  quite  close  to  any  point  on  the  unit 
circle.  Thus,  although  the  point  1  remains  fixed  under  iterates  of  f,  there 
always  exists  a  point  close  by  that  will  be  moved  somewhere  else  under  itera¬ 
tion.  A  similar  fate  falls  upon  all  points  in  {|z|  =  1},  and  thus  it  is 
possible  to  see  that  the  Julia  set  of  f(z)  =  z2  is  { | z J  =  1}. 

By  adding  a  constant  X  with  relatively  small  absolute  value,  the  dynamical 
system  produced  by  f(z)  =  z2  +  X  will  still  have  a  Julia  set  that  is  a  simple 
closed  curve.  However,  this  curve  exhibits  a  quasi-sel f-si mi lar i ty ,  i.e.,  it 
is  a  fractal.  As  the  value  of  jx|  gets  larger,  the  curve  degenerates  and  no 
longer  has  a  nicely  defined  inside  and  outside.  (Fig.  5  shows  various  Julia 
sets  produced  by  different  functions.) 

This  behavior  shows  up  in  the  parameter  space  image  of  +  \  and  is 
represented  by  the  Mandelbrot  set.  Recall  that  this  set  is  generated  by 
fixing  7  and  varying  X.  By  definition,  the  Mandelbrot  set  is  the  set  of 
complex  numbers  for  which  the  dynamical  system  generated  by  f(z)  =  z2  +  A  has 
a  connected  Julia  set.  For  X  =  0,  the  Julia  set  is  {|z|  =  1}.  As  |x|  in¬ 
creases,  the  resultant  Julia  set  of  the  dynamical  system  generated  by  f  will 
degenerate  from  a  simple  closed  curve.  When  |x|  gets  sufficiently  large,  the 
attractive  basis  inside  the  curve  bifurcates,  i.e.,  splits.  When  this  occurs, 
the  boundary  of  the  Mandlebrot  set  has  been  reached.  Theoretically,  it  has 
been  proven  that  this  behavior  is  completely  determined  by  the  growth  of  z  =  0 
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under  iteration  of  z2  +  X.  If  j  f n ( 0 ) |  is  sufficiently  large,  say  larger  than 
2,  then  fn(0)  will  converge  to  «.  This  indicates  that  the  Julia  set  for  that 
value  of  X  is  not  a  simple  closed  curve.  However,  if  J  f n  C  0 ) |  remains  bounded, 
then  the  Julia  set  is  a  simple  closed  curve. 

2. 2  Algor i thms 

The  programs  in  listing  1  allow  the  user  to  plot  some  pointillistic 
fractal  images  called  Julia  sets  (listings  are  provided  at  the  end  of  this 
section).  They  work  for  all  values  of  X  and  any  starting  point.  They  require 
only  a  simple  point  plotter  to  produce  images,  and  only  as  much  memory  as  the 
number  of  iterates  desired.  The  programs  were  written  in  VAX-11  FORTRAN  (DEC) 
because  it  supports  a  complex  variables  format  and  has  many  intrinsic  math 
functions.  If  FORTRAN  is  not  available,  the  code  listings  can  serve  as  a 
model  for  writing  a  program.* 

(For  the  reader  who  has  braved  the  theory  section,  an  explanation  of 
why  these  algorithms  work  is  now  simple.  First,  the  Julia  set  is  the  set  in 
which  the  function  does  not  have  convergent  neighborhoods.  It  is  preserved 
under  backward  iteration,  (f-1)n.  For  example,  if 

w  =  f(z)  =  z2  +  X  , 

z  =  f"'(w)  =  ±  /w  -  X  , 

and  so 

(f"')n  =  f- 1  0  ...  0  f~l  ,  n  times. 

*  Without  t',e  complex  variable  format ,  support  routines  have  to  be  written.  For  example,  multiplication 
must  follow  the  r  de 

ta,b)  *  ( c,d )  =  (ac  -  bd,  ad  +  be)  . 

A  comfdex  square  root  can  be  written  using  polar  coordinates.  Since 

V7  =  Vrc“  =  4r  , 

by  an  application  of  Euler' s  formula  and  the  half-angle  formulas,  we  obtain 

r  rT. J 1  +  cos  t\m  ,  .  (  !  -  cos  t\m 

*  Vr  ±(.  I  )  V  3  )  ' 


where  signs  are  chosen  according  to  quadrant.  Therefore,  because 

|  z  |  =  r  =  (x2  +  v2)112  ,  x  =  r  cos  t  ,  and  y  =  r  sin  t  , 

we  obtain 

l x , 


1  2 


Therefore,  to  get  an  iterate  into  the  Julia  set,  iterate  backwards,  choosing  a 
branch  (+)  of  the  square  root  at  random.  Thus,  the  iterate  is  never  allowed 
to  converge,  and  therefore  must  land  in  the  Julia  set.  Once  captured  there, 
it  just  moves  around  within  the  set.) 

Color  dynamical  images  in  C  are  produced  by  a  different  procedure 
(see  fig.  6,  center  spread,  pp  22  and  23).  Here,  the  colors  represent  conver¬ 
gence  rates.  After  a  domain  has  been  chosen,  the  image  is  produced  by  pro¬ 
ceeding  pixel  by  pixel  across  that  domain,  iterating  the  function  for  the 
value  of  z  represented  by  that  pixel.  Iteration  continues  until  j  f n ( z ) | 
reaches  a  certain  size,  or  the  function  has  been  iterated  a  predetermi ned 
maximum  number  of  times. 

The  number  of  iterations  then  determines  the  color  of  that  pixel. 
(This  is  essentially  the  same  procedure  which  was  outlined  by  Dewdney,  1985.) 
It  is  interesting  to  note  that  in  producing  these  color  images,  it  is  possible 
to  see  the  various  orbits  of  the  regions  computed. 

The  color  images  in  parameter  space  are  produced  in  the  same  way  as 
described  in  the  procedure  above.  However,  in  this  situation,  a  given  value 
of  z  is  fixed  throughout  the  entire  procedure,  while  the  position  of  the  pixel 
determines  the  value  of  A. 

The  color  images  of  both  the  C  dynamics  and  the  parameter  spaces  are 
not  limited  to  iterating  f(z)  =  z2  +  A.  All  the  color  images  shown  in  figure 
6  came  from  a  different  dynamic. 

Listing  2  is  a  FORTRAN  program  for  iterating  Newton’s  method  of 
finding  zeroes  as  applied  to  f(z)  =  z7  -  1  .  This  produced  the  pixel  informa¬ 
tion  to  create  figure  6(i). 

There  were  a  few  tricks  involved  in  computing  and  storing  the  data 
files  for  these  color  images.  As  this  type  of  program  is  computationally 
intensive  and  requires  quite  a  bit  of  computer  time,  computations  were  simpli¬ 
fied  when  possible  (e.g.,  using  |z|2  instead  of  |z|,  thus  eliminating  a  square 
root  for  each  iteration).  Since  the  resulting  data  file  would  occupy  much 
disk  space,  only  the  minimal  amount  of  information  needed  to  produce  an  image 
was  stored.  Pixels  were  computed  sequentially  on  a  given  row  with  the  func¬ 
tion  producing  a  color  value  for  each  pixel.  Adjacent  pixels  of  the  same 
color  were  considered  a  horizontal  vector.  When  a  new  color  value  for  a  pixel 
was  computed,  the  previous  pixel's  column  position  (terminal  point  for  the 
vector)  and  the  vector's  color  were  loaded  into  a  large  data  buffer: 

databuf(n)  «-  x 
databuf(n*1)  <-  oldcourit. 
n  ♦  n  ♦  2 
oldcount  ♦  count 

X  «-  X  ♦  1 

When  the  buffer  became  full,  it  was  written  to  a  disk  file  as  i  block 
of  unformatted  data.  A  second  program  i  listing  3)  must  be  used  to  plot  the 
vector  file.  Other  techniques  that  are  hardware  dependent  were  also  used. 


c  Plots  Julia  sets  for  quadratic  maps  by  iterating 
c  w  -  SQRT( z  -  c). 

c  Some  subroutine  calls  are  intrinsic  to  VAX-11  FORTRAN  (DEC) 
c  Programmer  :  S.  Casey 


PROGRAM  FRAC 
COMPLEX* 16  Z.  C.  CDSQRT 
REAL* 8  x,  y 

INTEGER *4  Niter,  I seed,  Tiaeseed 

WRITE  (6,*)  '  This  routine  plots  Julia  sets  for  f(z)  -  z**2  +  c 
WRITE  (6,*)  '  Enter  the  constant  ...  c-(a,b)' 

READ  ( 5 , * )  c 

WRITE  (6.*)  '  Enter  the  initial  z  value  ...  z-(a,b)' 

READ  ( 5 , * )  z 

WRITE  (6,*)  '  Enter  the  #  of  backward  iterates  ...  Niter' 

READ  (5.*)  Niter 

Iseed  -  Timeseed( )  !  get  seed  for  random  number  generator 

c  RAN( Iseed)  returns  a  floating-point  number  >-  0.0  and  <  1.0 

OPEN  (  10.  FILE- ' FRAC . DAT ' ,  STATUS- ' NEW ' ,  ERR-999,  IOSTAT-IOS, 

1  CARRIAGECONTROL* ' LIST '  ) 

DO  I  =  1 ,  Niter  !  Plots  positive  branch 
IF  (  RAN( Iseed)  .LT.  0.5  )  THEN 
z  -  -1 .  *  z 
ENDIF 
z  -  Z  -  c 
z  -  CDSQRT ( z ) 
x  -  DREAL(z) 
y  -  DIMAG ( z ) 

IF  (  I  .GE.  11  )  THEN  !  Let  fn(z)  converge  into  Julia  set 
WRITE ( 10 , * )  x,  y 
ENDIF 
ENDDO 

DO  I  -  1 ,  Niter  !  Other  branch 
IF  (  RAN( Iseed)  .LT.  0.5  )  THEN 
z  -  -1 .  *  z 
ENDIF 
z  -  z  -  c 

z  -  -1.  *  CDSQRT ( z ) 
x  -  DREAL(z) 
y  -  DIMAG(z) 

IF  (  I  .GE.  11  )  THEN  !  Let  fn(z)  converge  into  Julia  set 
WRITE( 10 , * )  x,  y 
ENDIF 
ENDDO 

CLOSE( 10) 

STOP 

999  WRITE  (6,*)  '  Error  opening  new  file  FRAC  DAT  -  IOS 
STOP 
END 


INTEGER *4  FUNCTION  Timeseed  () 

This  function  returns  a  large,  odd  integer  to  serve  as  an  initial 
seed  for  a  random  number  generator. 


c 

c 

c 

c 


Timeseed  -  INT(SECNDS(0.0))  !  get  number  of  seconds  since  midnight 
IF  (  MOD(Timeseed,2)  . EQ.  0  )  Timeseed  -  Timeseed  +  1  !  odd  value 

RETURN 
END 


Plots  Julia  sets  for  f(z)  -  z**2  +  L*z  by  iterating 
l/2(  -L  (+/-)  (L**2  +  4*z)  ). 

Some  subroutine  calls  are  intrinsic  to  VAX-11  FORTRAN  (DEC) 
Programmer  :  S.  Casey 


PROGRAM  FRAC2 

COMPLEX* 16  z,  L,  Det ,  Rootunity,  CDSQRT 
REAL* 8  x,  y 

INTEGER *4  K,  Niter,  Iseed,  Timeseed 
CHARACTER  Answer 

WRITE  (6,*)  '  This  routine  plots  Julia  sets  for  f(z)  -  z**2  +  L*z.' 

WRITE  (6,*)  '  The  results  are  interesting  if' 

WRITE  (6,*)  '  L  is  a  root  of  unity.' 

WRITE  (6,*)  '  Do  you  wish  to  enter  a  root  of  unity?' 

WRITE  (6.*)  '  Answer  y  for  yes. ' 

FORMAT  (A) 

READ  (5,95)  Answer 

IF  ((Answer  . EQ.  'Y')  .OR.  (Answer  . EQ.  'y'))  THEN 
WRITE  (6,*)  '  Enter  the  integer  denominator.' 

READ  (5, *)  X 
L  =  Rootunity  (X) 

WRITE  (6, *)  '  L  -  ' ,  L 
GO  TO  20 
ENDIF 

WRITE  (6,*)  '  Enter  the  constant  ...  L=(a,b)' 

READ  ( 5 , * )  L 

WRITE  (6,*)  '  Enter  the  initial  z  value  ...  z=(a,b)' 

READ  ( 5 , * )  z 

WRITE  (6,*)  '  Enter  the  #  of  backward  iterates  ...  Niter' 

READ  (5,*)  Niter 

Iseed  -  Timeseed( )  !  get  seed  for  random  number  generator 

RAN( Iseed)  returns  a  floating-point  number  >-  0.0  and  <  1.0 

OPEN  (  10.  FILE- ' FRAC2 . DAT ' ,  STATUS- ' NEW ' ,  ERR-999,  IOSTAT-IOS, 

'  CARRIAGECONTROL- ' LIST '  ) 


o  o 


Listing  1.  Code  producing  data  for  Julia  set  images  (cont'd) 


999 


C 


c 


DO  I  -  1,  Niter 

DET  -  CDSQRT  (  (L  *  L)  +  (4.  *  z)  )  | 

IF  (  RAN (I seed)  .LT.  0.5  )  THEN  ) 

DET  -  -1.  *  DET  I 

ENDIF  j 

z  -  (1./2.)  *  (  (-1.*  L)  +  DET  ) 

X  -  DREAL(z)  I 

y  -  DIMAG(z) 

IF  (  I  .GE.  11  )  THEN  !  Let  fn(z)  converge  into  Julia  set 

WRITE(10,*)  x.  y  ] 

ENDIF 

ENDDO  i 

CLOSEC 10)  ! 

STOP 

WRITE  (6,*)  ’  Error  opening  new  file  FRAC2.DAT  -  IOS  ' 

STOP 

END 


INTEGER *4  FUNCTION  Timeseed  () 

This  function  returns  a  large,  odd  integer  to  serve  as  an  Initial 
seed  for  a  random  number  generator. 

Timeseed  -  INT(SECNDS(0.0))  !  get  number  of  seconds  since  midnight 
IF  (  MOD (Time seed. 2)  .  EQ.  0  )  Timeseed  -  Timeseed  +  1  !  odd  value 

RETURN 
END 


COMPLEX* 16  FUNCTION  Rootunity(K) 

COMPLEX* 16  CDEXP ,  ITPIK 
REAL* 8  PI,  TPIK 
INTEGER *4  X 

PI  -  3.14159265358979323846 
TPIK  -  (2.  *  PI)  /  K 
ITPIK  -  (0.0, 1.0)  *  TPIK 
Rootunity  -  CDEXP(  ITPIK  ) 

RETURN 

END 
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Listing  2.  Code  producing  pixel  information  for  iteration  of  Newton's 
method  as  applied  to  f(z)  =  z7  -  1 

c - 

c  Iteration  of  Newton's  method  on  f(z)  -  z**7  -  1. 
c  Some  subroutine  calls  are  intrinsic  to  VAX-11  FORTRAN  (DEC) 
c  Programmers:  R.  Miller  and  S.  Casey 

c - 

PROGRAM  NEWT 

COMPLEX* 16  Znew,  Zold.  Ztemp 
REAL* 8  Diffabs 

INTEGER *2  iterations,  oldcount.  X,  count,  rows,  cols,  M,  N 
INTEGER  *  2  dat  abuf (16384) 

INTEGER *4  K,  J 

REAL* 8  acorner,  bcorner,  side,  gap,  realc,  image 
COMMON  / BLOCK 1/  dat abuf,  K 

c  Farms  passed  from  command  line:  iterations,  acorner,  bcorner,  side,  cols 

c - 

c  This  routine  calculates  pixel  information  for  a  rectangular  region  R. 
c  iterations  -  maximum  number  of  times  calculating  loop  is  executed 
c  acorner  -  Real  part  of  coordinate  of  lower  left  hand  corner  of  R 

c  bcorner  -  Imaginary  part  of  coordinate  of  lower  left  hand  corner  of  R 

c  side  -  length  of  horizontal  side  of  R 
c  rows, cols  -  number  of  pixels  in  R 

c - 

ACCEPT  *,  iterations,  acorner,  bcorner,  side,  cols 
rows  -  cols 

OPEN(  10,  FILE- ' NEWT . VEC ' ,  STATUS- ' NEW ' .  IOSTAT-IOS, 

*  ERR-999 ,  FORM- ' UNFORMATTED '  ) 


CALL  OUTBUF (  cols,  iterations  ) 

gap  -  side  /  REAL(rows) 

DO  N  -  rows  -  1,  0,  -1 

image  -  REAL(N)  *  gap  +  bcorner 
DO  M  -  0.  cols  -  1 

realc  -  REAL(M)  *  gap  +  acorner 
Znew  -  DCMPLX(  realc,  image  ) 
count  -  0 
Diffabs  -  1.0 


load  first  data  pair 


compute  pixels 


combine  real /imaginary  parts 


DO  WHILE ((count  .LT.  iterations ). AND . (Diffabs  .GT.  0.0001)) 
Zold  -  Znew 
Ztemp  -  7.0  *  Znew* *6 

IF  (  CDABS( Ztemp)  .LT.  .00001  )  THEN  !  absolute  value 
count  -  0 
GOTO  222 
ENDIF 

Znew  -  ((6.0  *  Znew* *7)  +  1.0)  /  Ztemp 
Diffabs  -  CDABS(  Znew  -  Zold  ) 
count  -  count  +  1 
ENDDO 
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Listing  2.  Code  producing  pixel  information  for  iteration  of  Newton's 
method  as  applied  to  f(z)  =  z7  -  1  (cont’d) 


Load  vector  data:  column  X  or  M  and  color  count'  or  oldcount 

52  IF(  M  .  EQ.  0  )  THEN 

oldcount  -  count 

X  -  0 

ELSE  IF(  M  . EQ.  (cols  -  1)  )  THEN 

IF(  count  .NE.  oldcount  )  CALL  OOTBUF(  X,  oldcount  ) 
CALL  OUTBUF (  M,  count  ) 

ELSE 

IF(  count  .HE.  oldcount  )  THEN 
CALL  OUTBUF (  X.  oldcount  ) 
oldcount  -  count 
ENDIF 
X  -  X  +  1 
ENDIF 

ENDDO 

ENDDO 

IF( (  K  .GT.  0  )  .AND.  (K  .LT.  16384))  THEN 
DO  J  -  K+l .  16384 

databuf (J)  -  0!  fill  remainder  of  buffer 
ENDDO 

WRITE(IO)  databuf!  write  last  record 
ENDIF 

CLOSE ( 10) 

CALL  EXIT 

19  WRITE  (6.*)  '  Error  opening  new  file  NEWT.VEC  IOS 
CALL  EXIT 
END 


SUBROUTINE  OUTBUF(  X.  color  ) 

INTEGER *2  X.  color 

INTEGER *2  databuf ( 16384 ) 

INTEGER*4  K/0/ 1  data  buffer  index 

COMMON  / BLOCK 1/  databuf,  K 

K  =  K  +  1 

databuf (K)  -  X 
K  -  K  +  1 

databuf(K)  -  color 
IF(  K  .  EQ.  16384  )  THEN 

WRITE(IO)  databuf  1  write  vector  data 

K  -  0 

ENDIF 


RETURN 


Pseudocode / FORTRAN  listing  illustrating  how  to  plot  vector  files. 

The  variable  'iterations'  can  be  used  for  loading  color  look-up  table 
Square  images  are  produced  (  rows  -  cols  ) . 

Programmer  :  R.  Miller 


INTEGER  *  2  iterations.  X.  Y.  color, 
INTEGER *2  xcenter,  ycenter 
CHARACTER  filename *40 

rows,  cols 

filename  -  'NEWT.VEC' 

OPEN(  10.  FILE-filename,  STATUS-  OLD' .  FORM- ' UNFORMATTED '  ) 

CALL  READBUF(cols ,  iterations)  ! 
rows  -  cols 

get  first  data  pair 

allocate  I/O  device 

enter  graphics  mode 
reset  graphics  device 

xcenter  -  cols  /  2 
ycenter  -  rows  /  2 
set  screen  and  coordinate  origins 

to  center  image  on  the  screen 

load  'iterations'  number  of  colors  into  look-up  tables 

Y  -  rows  -  1  1 

cols  -  cols  -  1  > 

range  of  rows  for  plot:  0  to  rows 
range  of  columns  :  0  to  cols 

plot  left  to  right  and  top  to  bottom  of  screen 

DO  WHILEC  Y  .GE.  0  ) 

MOVE  0 . Y  ! 
CALL  READBUF(X,  color)  ! 
DO  WHILE(  X  .LT.  cols  ) 

Move  to  beginning  of  row  Y 

get  first  vector  data  pair  for  row 

VALUE  color  ! 
DRAW  X+l.Y  ! 
CALL  READBUF ( X ,  color)  ! 

set  current  drawing  color 

draw  to  X+l  to  avoid  a  MOVE  X+l.Y 

get  next  vector  data  pair 

ENDDO 

VALUE  color  ! 

DRAW  X , Y  ! 

Y  -  Y  -  1 

set  drawing  color 

draw  last  vector  of  row  Y 

ENDDO 

exit  graphics  mode 
de-allocate  I/O  device 

CLOSE ( 10 )  ' 

close  file 

CALL  EXIT 

Listing  3-  Code/pseudocode  for  plotting  color  fractal  images  (cont'd) 


SUBROUTINE  READBUF(  X,  color  ) 

INTEGER *2  X,  color 
INTEGER *2  databuf( 16384) 

INTEGER*4  K/0/  !  data  buffer  index 

IF(  X  .  EQ.  0  )  READ(IO)  databuf 

X  -  X  +  1 
X  -  databuf(X) 

X  =  X  +  1 

color  -  databuf(X) 

IF(  X  .  EQ .  16384  )  X  =  0 

RETURN 

END 


3.  SEED  FRACTALS 

Fractal  images  may  also  be  generated  by  the  repetition  of  a  given  geomet¬ 
ric  pattern.  Examples  of  these  types  of  fractals  are  seen  in  figures  1  to  *1. 
These  images  are  generally  called  seed  fractals. 

By  one  definition,  a  fractal  curve  is  a  curve  for  which  the  fractional 
dimension  exceeds  the  topological  dimension.  Unless  the  curve  is  a  space¬ 
filling  Peano  curve,  this  topological  dimension  is  1 .  In  some  instances,  the 
curve  may  be  a  simple  closed  curve,  as  in  the  case  of  Koch's  snowflake.  If 
the  fractal  is  a  seed  fractal,  the  curve  is  constructed  by  a  limit  procedure 
where  a  given  seed  design  (some  geometric  shape)  is  scaled  and  repeated.  This 
produces  a  curve  that  is  self-similar;  that  is,  the  curve  pattern  repeats 
itself  on  any  level  of  magnification. 

3,1  Procedure  to  Produce  Seed  Fractal  Images 

All  the  curves  in  figures  1  to  4  are  seed  fractals.  In  each  of 
these,  the  algorithm  to  produce  them  was  a  twofold  process.  The  basic  pattern 
had  to  be  calculated,  scaled,  and  moved  to  its  proper  place  via  a  similarity 
transformation.  Simultaneously,  a  data  structure  had  to  be  set  up  in  order  to 
prepare  for  the  next_  level  of  iteration.  The  first  step  was  achieved  through 
trigonometry  and  linear  algebra.  The  second  step  was  handled  through 
counting. 


Listings  U  and  5  are  the  code  listings  for  Sierpinski's  gasket  (fig. 
7)  and  the  Cantor-Lebesgue  function  (fig.  2).  The  first  of  these  demonstrates 
the  drawing  process,  while  the  second  provides  an  example  of  the  counting 
process.  (Both  listings  are  given  at  the  end  of  this  section,  pp  28  to  31. 


3.2  Calculation  of  Topological  and  Fractional  Dimension 

Seed  fractals  present  a  good  opportunity  to  demonstrate  the  calcula¬ 
tion  of  fractional  dimension,  as  is  seen  in  the  following  examples.  Fraction¬ 
al  dimension  is  a  precise  gauge  on  how  much  an  object  "wiggles  about,"  and  is 
given  by  a  real  number.  It  is  different  from  our  more  intuitive  understanding 
of  dimension,  which  is  expressed  precisely  by  topological  dimension.  (Back¬ 
ground  for  this  section  is  provided  by  Guckenheimer  and  Holmes  (1983), 
Hurewicz  and  Mailman  (1948) ,  Lehto  and  Virtanen  (1973),  and  Munkres  (1973)-) 

Intuitively,  given  a  mathematical  object  in  Euclidean  3_space,  that 
object  is  usually  thought  of  as  having  dimension  0,  1 ,  2,  or  3  (which  are  the 
dimensions  of  a  point,  line  segment,  square,  and  cube,  respectively).  This 
intuitive  dimension  is  topological  dimension.  Topological  dimension  is  always 
given  by  an  integer,  and  corresponds  to  the  minimal  integer  value,  say  m,  for 
which  the  following  holds:  given  a  topological  space  X,  and  an  open  cover  A 
of  that  space,  there  is  a  refinement  B  of  A  that  has  order  m  +  1  .  Here,  a 
collection  B  of  subsets  of  A  has  order  m  +  1 ,  if  some  point  of  A  lies  in  m  +  1 
elements  of  B,  and  no  point  of  A  lies  in  more  than  m  +  1  elements  of  B. 

Let  the  diameter  of  a  set  be  the  least  upper  bound,  or  supremum  (sup) 
of  the  set  of  distances  between  points  in  the  set. 

Example:  The  unit  interval  I  =  [0,1]  has  topological  dimension  1. 

Let  I  be  endowed  with  the  subspace  topology  inherited  from  the  Eu¬ 
clidean  metric  topology  on  R. 

To  prove  that  dim  1*0,  assume  that  the  unit  interval  is  connected, 
that  is,  does  not  have  two  disjoint  nonempty  subsets  whose  union  equals  I. 
For  0  <  e  <  1,  let  A  be  any  open  covering  of  sets  of  diameter  less  than  e. 

Suppose  that  A  has  order  1.  Then,  no  two  elements  of  A  intersect.  Also, 

since  e  <  1 ,  A  must  contain  at  least  two  elements.  Let  U  be  one  element  of  A, 

and  let  V  be  the  union  of  the  others.  Then  U  fl  V  =  0  and  U  U  V  =  I,  contra¬ 

dicting  connectedness. 

Next,  let  A  be  any  open  cover  of  I.  Then, 

A  =  {[0,32)1  ...  ,  (an-2  > 1 J 1 

for  some  partition  of  the  unit  interval. 

0  —  3q  K  a  ^  K  ...  ^  an_ j  ^  —  1  • 

Let  c  =  min  (|ai  -  a^J}.  Now,  refine  P,  to  a  partition  P2  so  that 
P2  =  {bi:  0  =  bQ  <  b,  <  ...  <  bm-1  <  bm  =  1) 
with  | bA  -  bj_, |  <  e/2  for  all  i. 
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Figure  6.  Color  images  from  dynamics  in  C:  (a)  f(z)  =  X*exp(z),  parameter 
space;  (b)  f(z)  =  A*exp(z2),  parameter  space  (note  symmetry,  as  expected); 

(c)  f(z)  =  A-sin(z),  parameter  space  (note  "mi ni -Mandelbrot  sets"); 

(d)  f(z)  =  A*cos(z),  parameter  space;  ( e )  f(z)  =  A*tan(z),  parameter  space; 

(f)  f(z)  =  A'tan(z),  parameter  space  (enlargement  of  a  region  in  (e)); 

(g)  f(z)  =  z"  -  z  -  A,  parameter  space;  (h)  f(z)  =  z4  -  z  -  A,  parameter 
space  (enlargement  of  a  region  in  (g));  (i)  Newton’s  method  applied,  f(z)  = 
z7  -  1,  C-dynamical  system;  and  (j)  f(z)  =  z2  +  A,  parameter  space 
(enlargement  of  a  region  on  the  boundary  of  the  Mandelbrot  set,  containing  an 
outline  of  the  original  set). 

Pixel  colors  represent  convergence  rates.  In  all  the  figures  above,  this 
rate  was  an  "escape  rate"--how  quickly  the  iterates  of  that  pixel  converged 
to  infinity.  In  general,  the  color  schemes  went  according  to  color  frequen¬ 
cy.  Thus,  red  was  slow  convergence,  and  violet  was  fast  convergence.  How¬ 
ever,  this  scheme  was  not  strictly  adhered  to.  Also,  scaling  was  required, 
especially  in  figures  (a)  to  (f)  (expotential  growth  -»  logarithmic  scaling). 

All  the  color  images  seen  were  produced  by  R.  Miller,  who  worked  on  compu¬ 
ter  graphics,  and  S.  Casey,  who  did  mathematical  programming.  The  images 
were  produced  in  the  spring  and  summer  of  198b. 


7.  Sierpin, 
dure  followe' 
ing  procedur 


,>l,t 


,vt.l 


t 

♦  . 

*  f 

t-vJ . 

.*  ■ 

*  •  f 

‘  ./  , • 

*  '  V  ! 

*  •/  *  * 

•V  '/J'. 

*-v  - 

,Vv-  h 

L- _ hi  a-  » 

-1.00  -  0.75  -0.50  -  0.25  0  00  0  25  0  50  0.75  1.00 

X-axis 


-100  -  0  75  -0  50  -0.25  0  00  0  25  0  50  0  75  1.00 

X-axis 


B  =  |  C 0 , b 2 ) .  ^3) .  •  •  .  .  ( bm-2  > 1  ]  J 

is  a  refinement  of  A  of  order  2. 
Since  this  is  true  for  any  cover  A, 
it  follows  that 

dim  I  £  1 

Therefore,  since  topological  dimen¬ 
sion  is  integer  valued,  dim  1=1.  O 

To  prove  that  the  topologi¬ 
cal  dimension  of  a  simple  (noninter¬ 
secting)  continuous  curve  is  equal  to 
1  requires  a  bit  more  machinery. 

In  order  to  calculate  frac¬ 
tional  dimension,  the  idea  of  cover¬ 
ing  a  set  by  open  sets  is  used.  How¬ 
ever,  in  this  situation,  the  diameter 
of  the  sets  used  plays  an  important 
role.  The  idea  of  fractional  dimen¬ 
sion  goes  back  to  the  nineteenth  cen¬ 
tury,  to  F.  Hausdorff  and  I.  Besoco- 
vitch.  The  following  is  the  formal 
definition  of  the  frs’tional  dimen¬ 
sion  named  after  them. 

Let  X  be  some  set  in  Euclid¬ 
ean  3-space.  Consider  all  coverings 
of  X  by  countably  many  sets  of  diam¬ 
eter  less  than  some  d  >  0.  Given  one 
such  covering,  say  [o^  ,  a,  . ..},  as¬ 
sociate  with  that  cover  the  number 
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Figure  7  (cont'd).  The  arrowhead: 

(g)  defining  procedure  followed  three 
times,  (h)  four  times,  and  (i)  five 
times . 


where  dp  is  the  diameter  of  the  set 
an  and  f>  >  0.  For  every  6  let 


i*(X)  =  lim  [inf  (£  d  ]] 


where  inf  (infimum)  is  the  greatest 
lower  bound.  The  quantity  u*  is 
called  the  B-dimensional  Hausdorff- 


«  *"•  -  '  -  -  *  -/  *.*  "S.  A  A,**  .  • 


Be?  covitch  outer  measure.  Then,  the  Hausdorff-Besocovitch  dimension  of  X  is 
given  by 


dim  X  =  inf  { B I u*(X)  =  0}  . 

p 

Again,  this  definition  is  best  seen  in  the  following  examples. 
Example — the  Cantor  middle  thirds  set:  Let  Cq  =  [0,  1].  Then, 
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To  get  an  upper  bound  on  dim  C, 
can  be  covered  by  2n  sets  of  length  3~n. 


first  work  with  the  set  Cn> 
Now, 


This  set 


and  so  (nonrigorously ) 


n 


-  2n(3'n)e  , 


u*(C)  =  lim  inf  2n(3~n)8 
8  n+» 


To  calculate  the  dimension  of  C,  one  must  find  the  inf  {bIuJ(C)  =  o}.  As 

B 

2n3  n8  =  x  «*  n(log  2  -  B  log  3)  =  log  x  , 


and 


log  x  +  ®  as 


n  ->  00 


B  < 


log  2 
log  3 


we  obtain 


dim  C  < 


log  2 
log  3  ’ 


With  a  bit  more  work ,  one  can  see  that  this  bound  is  sharp.  O 
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Example-  Koch*s  snowflake:  At  the  ri  stage  of  its  construction ,  the  snow- 
flake  can  be  covered  by  sections  of  length  3"n.  Thus, 


Therefore,  since 


we  obtain 


u*(K)  =  lim  inf  *in3  nB 


4ne  nB  =  x  «+  n(log  4  -  B  log  3)  =  log  x  , 


log  x  <  ®  as  n  «  8  < 


log  4 
log  3  , 


dim  K  <  -1  .  O 

-  log  3 


Discussions  of  topological  dimension  may  be  found  in  Munkres  (1975) 
and  Hurewicz  and  Wallman  (1998).  Discussions  of  fractional  dimension  may  be 
found  in  Lehto  and  Virtanen  ( 1 97  j )  and  Mandelbrot  (1983).  One  of  MandelDrot's 
main  underlying  themes  is  that  a  fractal  is  character i zed  by  its  fractional 
dimension.  Thus,  numerous  calculated  examples  can  be  found  throughout  the 


Listing  4.  Code  producing  vector  data  to  draw  Sierpinski’s  gasket 

c - 

c  This  routine  draws  Sierpinski's  Gasket. 

c  Some  subroutine  calls  are  intrinsic  to  VAX-11  FORTRAN  (DEC) 
c  Programmer  :  S.  Casey 

c - 

PROGRAM  SGASKET 
REAL’S  Con,  Scale 

REAL’S  Xnode(7 , 730) ,  Ynode(7,730)  !  Centers  of  the  triangles... 

REAL *8  Xreturn(3),  Yreturn(3) 

INTEGER  Iter,  I,  J,  X 

CHARACTER  *  1  Lntyp  !  Draws  either  points(P)  or  vectorsC V) . . . 

WRITE  (6,*)  '  Enter  the  #  of  iterates  ...  a  maximum  of  7.' 

READ  (5, *)  Iter 

OPEN  (10,  FILE= 'SGASKET.DAT' .  STATUS= ' NEW ' ,  ERR=999,  IOSTAT=IOS, 
’  CARRIAGECONTROL= ' LIST '  ) 

Con  =  3.0 

Con  -  DSQRT(Con) 

100  FORMAT( IX , F10 . 3 , IX , F10 . 3 , IX, A) 

Lntyp  =  ' P ' 

WRITE(10, 100)  0.0,  1.0,  Lntyp  !  Scaling  factors. 

WRITE(10, 100)  0.0,  0.0,  Lntyp 
Lntyp  -  'V' 

WRITE( 10. 100)  0.0,  0.0,  Lntyp  !  The  outside  triangle. 

WRITE(IO.IOO)  1.0,  0.0,  Lntyp 

WRITE(IO.IOO)  0.5,  (Con  /  2.0),  Lntyp 

WRITE(10, 100)  0.0,  0.0,  Lntyp 

Lntyp  =  'P'  !  Lifting  the  pen. 

WRITE(10. 100)  0.25,  (Con  /  4.0),  Lntyp 
Lntyp  =  'V' 

WRITE( 10, 100)  0.75,  (Con  /  4.0),  Lntyp  !  The  inside  triangle. 
WRITE( 10. 100)  0.5,  0.0,  Lntyp 
WRITEdO.  100)  0.25,  (Con  /  4.0),  Lntyp 

Xnode(l.l)  »  0.5 
Ynode(l.l)  -  Con  /  8.0 

DO  I  -  1,  Iter 

Scale  -  1.0  /  (2.0”(I+1)) 

DO  J  =  1,  3” ( 1-1 ) 

CALL  Gasketseed( Xnode( I , J) .  Ynode(I.J),  Scale, 

•  Xreturn,  Yreturn) 

DO  K  =  1.  3 

Xnode( ( 1  + 1 ) ,  (3*J-(3-K)))  «  Xreturn(K) 

Ynode( ( 1+ 1 ) ,  (3*J-(3-K)))  =  Yreturn(K) 

ENDDO 

ENDDO 

ENDDO 

CLOSE ( 10) 

STOP 

999  WRITE  (6.*)  '  Error  opening  new  file  SGASKET.DAT  -  ',  IOS 
STOP 
END 
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Listing  J4.  Code  producing  vector  data  to  draw  Sierpinski's  gasket  (cont'd) 


SUBROUTINE  Gasketseed(Xnode ,  Ynode,  Scale,  Xreturn,  Yreturn) 


REAL'S  Xnode,  Ynode,  Scale,  Con 
REAL *8  Xreturn(3),  Yreturn(3) 

REAL *8  Xvar(3,4) ,  Yvar(3,4) 

CHARACTER  *  1  Lntyp  !  Draws  either  points(P)  or  vectors(V).. 

100  FORMAT( 1X.F10. 3, IX , F10 . 3 , IX , A) 

Con  -  3.0 

Con  -  DSQRT(Con) 

Xreturn(l)  -  1.0  /  2.0 
Yreturn(l)  -  -Con  /  8.0 

Xreturn(2)  -  0.0 

Yreturn(2)  -  (3.0  *  Con)  /  8.0 


Xreturn(3) 

Yreturn(3) 


-1.0  /  2.0 
-Con  /  8.0 


DO  1=1,3 

Xreturn(I)  -  (2.0  *  Scale  *  Xreturn(I))  +  Xnode 
Yreturn(I)  =  (2.0  *  Scale  *  Yreturn(I))  +  Ynode 
ENDDO 


Xvar(I.l)  =  -1.0  /  2.0 
Yvar(I.l)  -  Con  /  4.0 

Xvar (1,2)  -  1.0  /  2.0 
Yvar(I,2)  -  Con  /  4.0 

Xvar (1,3)  >0.0 
Yvar(I,3)  -  -Con  /  4.0 

Xvar(I,4)  *  Xvar(I,l) 

Yvar(I,4)  =  Yvar(I,l) 

Lntyp  >  ' P ' 

DO  J  -  1 ,  4 

Xvar(I.J)  =  (Scale  *  Xvar(I,J))  +  Xreturn(I) 
Yvar(I.J)  =  (Scale  *  Yvar(I.J))  +  Yreturn(I) 
IF  (  J  .EQ.  2  )  Lntyp  -  'V' 

WRITE( 10 , 100)  Xvar(I,J),  Yvar(I.J),  Lntyp 
ENDDO 


ENDDO 


RETURN 


lv. 


l  4 1.  ilk'itk  ll.’lK  4I.  i|.  it.  it.  it.  it.  it^'iL ‘il. 'it.'il.'il. ‘it.'il.'il.’il.'il,1  |tt_  *1- *Al.  lln'tL ‘jt.1 
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Listing  5.  Code  producing  vector  data  to  draw  Cantor-Lebesgue  function 

c - 

c  This  routine  draws  the  Cantor-Lebesgue  funtion. 
c  Some  subroutine  calls  are  intrinsic  to  VAX- 11  FORTRAN  (DEC) 
c  Programmer  :  S .  Casey 


PROGRAM  CANTOR 

REAL’S  Xl( 8 , 256 ) ,  Yl(8,256),  X2(8.256).  Y2(8,256),  Xstep,  Ystep 

REAL’S  Xreturn/ 1 . 0/ ,  Yreturn/0 . 0/ .  Xscalel /0 . 0/ ,  Yscalel/0.0/ 

REAL’S  Xscale2/ 1.0/,  Yscale2/1.0/ 

INTEGER *4  Iter,  Count /0/ 

WRITE  (6,*)  '  Enter  the  #  of  iterates  ...  a  maximum  of  8. ' 

READ  (5, *)  Iter 

OPEN  (  10,  FILE= ' CANTOR . DAT ' ,  STATUS= ' NEW ' .  ERR=999,  IOSTAT=IOS 
CARRI AGECONTROL=  LIST'  ) 

DO  I  =  1,  Iter 

IF  (I  . EQ.  Iter)  THEN 

WRITE( 10 , * )  Xscalel,  Yscalel 
END  IF 

Count  =  Count  +  2”(I-1) 

DO  J  =  1 ,  Count 

Xstep  =  (1.0  /  (3. 0**1)) 

Ystep  -  (1.0  /  (2. 0**1)) 

IF  (J  . EQ.  1)  THEN 
X1(I+1,J)  =  Xstep 
YKI  +  1.J)  =  Ystep 
X2(  1+ 1 ,  J )  =  XKI+1.J)  +  Xstep 
Y2( 1+1 , J)  =  Y1(I+1,J) 

IF  (I  .EQ.  Iter)  THEN 

WRITE(10,*)  X1(I+1,J),  Y1(I+1,J) 

WRITEdO,*)  X2(  1  + 1 ,  J)  ,  Y2(I+1,J) 

END  IF 
END  IF 


IF  (MOD( J , 2) 
XKI  +  1 ,  J) 
YKI+I ,  J) 
X2(I+1 , J) 
Y2( 1+ 1 , J) 
IF  (I  . EQ 
WRITE( 
WRITE( 
END  IF 
ELSE 

XKI+1.J) 
YKI+I,  J) 
X2(I+ 1 , J ) 
Y2( 1+ 1 , J) 
IF  (I  . EQ 
WRITE( 
WRITE( 
END  IF 
END  IF 
ENDDO 


. EQ.  0)  THEN 
=  X1(I, (J/2)) 

-  Y1(I , ( J/2)  ) 

=  X2(I, (J/2)) 

=  Y2( I , ( J/2 )  ) 

.  Iter)  THEN 

10,  *)  XKI+1.J),  YKI+1.J) 
10. * )  X2( 1+ 1 , J) ,  Y2(I+1 , J) 


=  X2( I , ( J/2) )  +  Xstep 
=  Y2( I , (J/2))  +  Ystep 
=  X1(I+1,J)  +  Xstep 
=  Yld+l, J) 

.  Iter)  THEN 

10,  ’  )  XKI+1,  J).  Y1(I+1,J) 
10,*)  X2( 1+ 1 , J ) ,  Y2( 1+ 1 , J) 


Listing  5.  Code  producing  vector  data  to  draw  Cantor-Lebesgue 

function  (cont'd) 

IF  (I  . EQ.  Iter)  THEN 

WRITE( 10 , * )  Xscale2.  Yscale2 
WRITE(10,*)  Xreturn,  Yreturn 
END  IF 
BNDDO 

CLOSE(IO) 

STOP 

999  WRITE  (6,*)  '  Error  opening  new  file  CANTOR . DAT  -  IOS 
STOP 
END 

4.  SQUIG  FRACTALS 

As  is  known,  fractal  images  seem  to  imitate  nature  rather  well.  Certain 
types  of  fractals  employ  a  controlled  random  motion  to  produce  images  of 
landscapes,  trees,  etc.  Mandelbrot  calls  these  squig  fractals. 

Fractals  have  now  worked  their  way  into  computer  graphics  art,  and  even 
into  the  movie  theatre.  The  sharp  landscape  images  now  seen  in  some  films 

were  generated  by  taming  Brownian  (or  random)  motion.  The  techniques  for 

generating  these  squig  fractals  are  somewhat  similar  to  those  employed  in 
generating  the  seed  images.  Again,  an  image  must  be  drawn,  and  a  data  struc¬ 
ture  which  anticipates  the  next  level  of  iteration  must  be  in  place.  However, 
the  big  difference  is  that  in  the  production  of  squigs,  a  computerized  "coin 
flip,"  provided  by  calling  a  random  number  routine,  determines  the  final  shape 
of  the  fractal. 

Uncontrolled  random  motion  appears  too  rough  to  model  nature.  This  can  be 
seen  easily  by  randomly  choosing  (x,  y)  coordinates,  and  connecting  the  dots 

(see  fig.  8).  However,  if  some  control  is  introduced  into  this  process,  such 

as  scaling,  transforming,  or  filtering,  a  pattern  which  imitates  nature’s  own 
"controlled  randomness"  seems  to  emerge.  This  technique  of  creating  images  by 
controlled  random  motion  has  been  employed  successfully  in  the  generation  of 
landscape  images,  including  mountain  ranges  and  lakes  (see  fig.  9). 

Unfortunately,  the  computer  code  to  generate  even  the  most  basic  of  the 
squig  images  is  long  and  tedious  (listing  6,  pp  3^-40 ) . 

We  can  demonstrate  the  fundamental  idea  behind  the  creation  of  a  squig 
fractal  by  discussing  a  planar  squig  curve.  Given  a  simple  closed  polygonal 
curve,  the  squig  procedure  can  be  applied  to  produce  a  controlled  Brownian 
path  which  circumvents  nearly  the  same  area  as  the  initial  curve. 

Cover  the  initial  curve  with  a  grid  of  sufficiently  small  mesh  size  (say, 
d/9,  where  d  =  minimal  distance  between  vertices).  Let  the  length  of  one  line 
segment  in  the  grid  be  denoted  by  L.  Align  the  mesh  so  that  none  of  the 
polygon's  vertices  matches  a  vertex  of  the  mesh.  Thus,  in  each  square  that 


Figure  8.  Brownian  motion. 


X-axis 

Figure  9.  Skeletal  fractal  mountain. 


the  original  path  intersects,  the  polygonal  path  will  come  in  one  side  and 
exit  on  another.  Divide  the  box  into  four  equal  boxes  having  sides  of  length 
L/2.  At  this  point,  controlled  random  motion  will  produce  the  first  genera¬ 
tion  of  a  squig  curve. 

Consider  a  single  square.  Assume  that  the  path  enters  the  box  from  the 
left.  In  the  scale  of  the  finer  mesh  (L /<?'),  the  path  could  enter  the  box  from 
the  top  or  bottom.  In  the  squig  scheme,  this  is  determined  by  a  weighted  coin 
flip.  e.g. , 

ran  (  )  <  0.2b  =  3 / U  weight 


If  this  is  the  first  box  processed  on  this  level  of  iteration,  then  i  fair 
coin  flip  (i.e.,  ran  (  )  £  0.8)  will  determine  entry.  If  not,  then  entry 
position  is  inherited  from  the  previous  pox.  Similarly,  in  the  scale  of  the 
finer  mesh,  there  are  two  choices  of  exit  f mm  each  side.  This  is  decided  by 
a  weighted  coin  flip.  (This  exit  position  then  determines  entry  of  the  adja¬ 
cent  box.)  Between  entry  and  exit,  there  are  22  possible  paths  on  the  level 
of  the  finer  mesh.  There  are  eight  exiting  from  the  opposite  side,  arid  seven 
each  from  the  adjacent  sides.  The  choice  of  any  of  these  paths  comes  from  a 
couple  of  weighted  coin  flips,  whi  -h  ict-e-mi  ne  whether  or  r  -  *  t  he  path  turns 
or  goes  straight  isee  fig.  in). 


Throughout  the  process,  t  tie  coin  fl 
straighter  path.  The  higher  the  weight, 
lower  the  fractional  dimension  of  the  curve 


weighted  in  f  tvor  of  the 
straighter  the  path,  md  the 


Once  this  process  of  choosing  t he  path  on  the  first  level  of  iteration  has 
been  completed,  f  tie  curve  must  be  covered  ty  boxer,  w : *■  h  side..  rd  length  LM . 
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Figure  10.  Possible  two-dimensional 
squig  paths. 


in  turn,  each  of  these  must  be 
divided  in  equal  fourths,  and  the 
next  level  of  iteration  in  the  squig 
procedure  must  be  done.  Iteration 
continues  a  predetermined  number  of 
times . 

The  squig  procedure  is  by  no 
means  limited  to  curves.  The  idea 
of  using  chance  to  produce  fractal 
benavior  is  a  particularly  powerful 
one.  This  is  witnessed  in  The  Frac¬ 
tal  Geometry  of  Nature,  for  Mandel¬ 
brot  devotes  nearly  the  entire 
second  half  of  the  book  to  this  sub¬ 
ject.  One  particularly  interesting 
use  of  the  squig  procedure  is  in  the 
production  of  semi  random  tree  paths 
called  "graftals"  (see  fig.  11). 
Estvanik's  article  on  this  subject 
(1986)  is  particularly  useful. 


Producing  a  squig  is  only  the  first  step  towards  producing  a  realistic 
landscape.  Smoothing,  coloring,  shading,  shadowing,  and  numerous  other  more 
standard  techniques  from  computer  graphics  must  then  be  applied  to  produce  the 
finished  product,  a  realistic  image.  It  is  interesting  to  note,  however,  that 
because  the  squig  (and  other)  fractal  methods  are  so  powerful,  these  tech¬ 
niques,  once  considered  esoteric,  are  rapidly  being  accepted  into  the  main¬ 
stream  of  computer  graphics. 
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Listing  6.  Code  for  plotting  skeletal  fractal  mountains 


c  - 

c  This  program  draws  the  skeleton  of  a  fractal  mountain 
c  range  by  means  of  the  addition  of  random  weights  to 
c  the  y  values  of  a  triangular  lattice.  The  con- 

c  struction  of  this  lattice  goes  from  0  to  Iter, 

c  Programmer  :  S.  Casey 

c  - 

c 


PROGRAM  MOUNTAIN 

REAL* 8  X( 0 : 6 ,  24768).  Y(0:6,  24768)  !  Array  of  endpoints 

REAL *8  Mx( 24768 ) .  My(24768)  <  Array  of  midpoints 

REAL* 8  Randbound  !  The  degree  of  roughness 

REAL* 8  Randshft  !  Function  giving  random  number  with  random  sign 
INTEGER *4  I,  J,  Iter,  Level 

INTEGER*4  Depth(0:6),  Numpts(0:6),  Numlns(0:6).  Numtrg(0:6) 

WRITE  (6.*)  '  Enter  the  #  of  iterates.  Cannot  do  more  than  6.' 

READ  (5, *  )  Iter 

WRITE  (6,*)  '  Enter  the  degree  of  roughness.  This  is  a  number' 

WRITE  (6.*)  '  between  0  and  1.' 

READ  ( 5 , *  )  Randbound 

WRITE  (6.*)  '  Enter  the  initial  triangle.' 

DO  I  -  1 ,  3 

READ  (5, *  )  X(0.  I) .  Y(  0 ,  I) 

ENDDO 

Depth(O)  -  2 
Numpts(O)  -  3 
Numlns(O)  =  3 
Numtrg(O)  =  1 
DO  I  =  1,  Iter 

Depth(I)  =  Depth(I  1)  +  2**(I-1) 

Numpts(I)  =  (  Depth(I)  *  (Depth(I)  +  1)  )  /  2 
Numlns ( I )  =  3  •  (  (  (Depth(I)  -  1)  *  Depth(I)  )  /  2  ) 

NumtrgC I )  -  4’ *  I 
ENDDO 

Level  -  0 
DO  I  ’  1,  Iter 

Level  =  Level  *  1 

CALL  Midpointsf  X.  Y,  Randbound.  Level.  Mx,  My  ) 

CALL  Reassignt  X.  Y.  Level.  Mx ,  My  ) 

ENDDO 

OPEN  (  10.  FILE  MOUNTAIN  DAT  .  STATUS' ’ NEW ' ,  ERR'996. 

*  IOSTAT  IOS.  CARRIAGECONTROL  LIST  ) 

CALL  Draw  X  Y .  Iter  ) 

CLOSE!  10  ' 

STOP 

998  WRITE  (6.*)  Error  opening  new  file  MOUNTAIN . DAT  =  IOS 
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Listing  6.  Code  for  plotting  skeletal  fractal  mountains  (eont'd) 


c - 

c  --  This  function  returns  a  bounded  random  number 
c  -Randbound/2* ‘Level  <=  num  <=  Randbound/2“Level 
c  with  random  sign, 

c 


REAL* 8  FUNCTION  Randshft!  Randbound,  Level  ) 

REAL* 8  Randbound  !  Variable  received  from  main 

REAL *8  Rand,  Sgn,  Scale  !  Local  variables 

INTEGER *4  Iseed,  Timeseed  !  For  the  random  number  generator 

INTEGER*4  Level  !  Level  of  iteration  in  main  routine 

LOGICAL  First/ . TRUE . /  !  Flag  which  insures  only  one  call  of  Timeseed 

IF  (  First  )  THEN 

Iseed  =  Timeseed! ) 

First  -  .FALSE. 

END  IF 

IF  (  RAN( Iseed)  .LT.  0.5  )  THEN 
Sgn  =  -1.0 
ELSE 

Sgn  =1.0 
ENDIF 


> 

V 

V 
■ 
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fi 
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Rand  =  RAN (I seed)  *  Randbound 
Scale  -  1.0  /  (2.0**Level) 

Randshft  =  Sgn  *  Scale  *  Rand 

RETURN 

END 


C - 

c  --  This  function  is  a  sleazy  trick  which  enables  people  with  typing 
c  disabilities,  like  me,  to  get  a  large  odd  integer  for  RAN( ) . 

c 

INTEGER *4  FUNCTION  Timeseed! ) 

Timeseed  =  INT! SECNDS! 0 . 0) ) 

IF  (  MOD! Timeseed , 2 )  . E Q.  0  )  Timeseed  -  Timeseed  +  1 

RETURN 

END 


C - 

c  --  This  subroutine  calculates  the  midpoints  of  the  lattice! level- 1 ) . 
c 

SUBROUTINE  Midpoints!  X.  Y.  Randbound,  Level.  Mi,  My  ) 


■  * 
I 


REAL'S  X!0:6,  24768),  Y!0:6,  24768) 

REAL'S  Randbound 

REAL'S  Mx! 24768).  My! 24768) 

REAL* 8  Randshft 

REAL *8  Rand.  XI,  X2,  Yl.  Y2 


Endpoints,  from  main 
Variable  received  from  main 
Returned  to  main 
Random  number  with  random  sign 
Local  variables 


v. 

w. 
w, 
>, 

V 
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Listing  6.  Code  for  plotting  skeletal  fractal  mountains  (cont'd) 


INTEGER *4  Level  !  Level  of  Iteration  in  main 

INTEGER *4  I,  J,  Depth,  Int ,  Count  !  Local  variables 
INTEGER *4  Leadpt(129).  Jump(0:129) 

Depth  -  2 

DO  I  -  1,  (Level-1) 

Depth  -  Depth  +  2**(I-1) 

ENDDO 


Leadpt ( 1 )  =1 
Int  -  0 

DO  I  =  2,  (Depth+1 ) 

Int  -  Int  +  1 

Leadpt(I)  -  Leadpt (1-1)  +  Int 
ENDDO 

Calculating  midpoints  of  line  segments  "slanting  left" 


Count  -  0 

DO  I  -  1.  (Depth-1) 

Jump(O)  -  0 

DO  J  -  1,  (Depth-I ) 

Jump(J)  -  (<7  +  I-  l)  +  Jump(J-l) 

XI  -  X(  (Level-1).  (  (Leadpt(I+l)-l)  + 

Y1  -  Y(  (Level-1),  (  ( Leadpt (I+l)-l)  + 

X2  -  X(  (Level-1),  (  (Leadpt(I+l)-l)  + 

Y2  -  Y(  (Level-1),  (  (Leadpt ( 1+1 ) -1 )  + 

Rand  -  Randshft(  Randbound,  Level  ) 
Count  -  Count  +  1 
Mx( Count)  -  (XI  +  X2)  /  2.0 
My( Count)  -  ((Yl  +  Y2)  /  2.0)  +  Rand 
ENDDO 
ENDDO 


Jump(J-l)  )  ) 
Jump(J-l)  )  ) 
Jump(J)  )  ) 
Jump(J)  )  ) 


Calculating  midpoints  of  line  segments  "slanting  right" 

DO  I  -  l.  (Depth-1) 

Jump(O)  -  0 

DO  J  -  1,  (Depth-I) 

Jump(J)  -  (J  +  I)  +  Jump(J-l) 

XI  -  X(  (Level-1).  ( Leadpt ( I )  +  Jump(J-l))  ) 

Yl  -  Y(  (Level-1),  (Leadpt(I)  +  Jump(J-l))  ) 

X2  -  X(  (Level-1),  (Leadpt(I)  +  Jump(J))  ) 

Y2  -  Y(  (Level-1).  ( Leadpt ( I )  +  Jump(J))  ) 

Rand  -  Randshft(  Randbound,  Level  ) 

Count  -  Count  +  1 
Mx( Count )  -  (XI  +  X2)  /  2.0 
My (Count )  -  ((Yl  +  Y2)  /  2.0)  +  Rand 
ENDDO 
ENDDO 


■ 
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Listing  6.  Code  for  plotting  skeletal  fractal  mountains  (cont'd) 

Calculating  midpoints  of  “horizontal"  line  segments 

DO  I  -  2.  Depth 

DO  J  -  Leadpt(I).  Leadpt(I+l )-2 
XI  -  X(  (Level-1).  J  ) 

Y1  -  Y(  (Level-l),  J  ) 

12  -  X(  (Level-1).  J+l  ) 

Y2  -  Y(  (Level-1),  J+l  ) 

Rand  -  Randshft(  Randbound,  Level  ) 

Count  -  Count  +  1 
to (Count )  -  (XI  +  X2)  /  2.0 
My (Count )  -  ((Yl  +  Y2)  /  2.0)  +  Rand 
ENDDO 

ENDDO 

RETURN 

END 

This  subroutine  reassigns  points  of  latticeC level-1 )  and  new 
midpoints  to  the  lattiee( level ) . 

SUBROUTINE  Reassign(  X,  Y,  Level,  to.  My  ) 

REAL  *  8  X(0:6,  24768),  Y(0:6,  24768)  !  Array  of  endpoints 

REAL *8  Mx( 24768) ,  My(24768)  !  Array  of  midpoints 


INTEGER *4  Level 

INTEGER *4  I,  J.  Int .  Jump,  Count 

INTEGER *4  Skipl,  Skip2 

INTEGER *4  Depth(0:7),  Leadpt(129) 


!  Level  of  iteration  in  main 
!  Local  variables 


Depth(O)  «  2 
DO  I  =  1 ,  Level +1 

Depth(I)  -  Depth(I-l)  +  2**(I-1) 

ENDDO 

Leadpt(l)  -  1 
Int  -  0 

DO  I  -  2,  Depth(Level+l ) 

Int  -  Int  +  1 

Leadpt(I)  -  Leadpt(I-l)  +  Int 
ENDDO 

Reassigning  old  lattice  points 

X(  Level.  1  )  -  X(  (Level-1),  1  ) 

Y(  Level.  1  )  =  Y(  (Level-1).  1  ) 

Jump  -  0 

DO  I  -  2,  Depth( Level- 1 ) 

Jump  -  Jump  +  1 

X(  Level.  Leadpt( I+Jump)  )  -  X(  (Level-1), 
Y(  Level,  Leadpt(I+Jump)  )  -  Y(  (Level-1), 
Skipl  -  0 

DO  J  -  Leadpt(I)+l ,  Leadpt ( 1+ 1 )- 1 
Skipl  -  Skipl  +  1 

X(  Level,  (Leadpt (I  +  Jump)  +  (2 * Skipl ) )  ) 
Y(  Level,  (Leadpt(I+Jump)+(2*Skipl ) )  ) 
ENDDO 
ENDDO 


Leadpt(I)  ) 
Leadpt(I)  ) 


X(  (Level-1), 
Y(  (Level-1), 


tow 
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Listing  6.  Code  for  plotting  skeletal  fractal  mountains  (cont'd) 

Reassigning  new  midpoints  to  lattice 

Reassigning  midpoints  of  line  segments  slanting  left 

Count  -  0 
Skipl  =  0 

DO  I  ■  1,  ( Depth( Level ) - 1 ) 

Skipl  -  Skipl  +  1 
Jump  -  0 

IF  (  MOD( Skipl,  2)  .NE.  0  )  THEN 
DO  J  -  1,  (Depth(Level)-I ) 

Jump  -(J+I-l)+  Jump 
IF  (  MOD( J .  2)  .NE.  0  )  THEN 
Count  -  Count  +  l 

X(  Level,  (Leadpt(I+l)-l+Jump)  )  -  Mx(  Count  ) 

Y(  Level,  (Leadpt(I+l)-l+Jumn)  )  -  My(  Count  ) 
ENDIF 

ENDDO 

ENDIF 

ENDDO 

Reassigning  midpoints  of  line  segemnts  slanting  right 
Skipl  -  0 

DO  I  -  1.  (Depth(Level)-l) 

Skipl  -  Skipl  +  1 
Jump  -  0 

IF  (  MOD( Skipl,  2)  .NE.  0  )  THEN 
DO  J  -  1.  (Depth(Level)-I) 

Jump  -  (J  +  I)  +  Jump 
IF  (  MOD( J,  2)  .NE.  0  )  THEN 
Count  «  Count  +  1 

X(  Level,  ( Leadpt ( I ) + Jump )  )  -  Mx(  Count  ) 

Y(  Level,  (Leadpt ( I )+Jump)  )  -  My(  Count  ) 

ENDIF 

ENDDO 

ENDIF 

ENDDO 

Reassigning  midpoints  of  horizontal  line  segments 

Skipl  -  1 
Skip2  -  0 

DO  I  =  2,  Depth( Level) 

Skipl  -  Skipl  +  1 

IF  (  MOD( Skipl,  2)  .NE.  0  )  THEN 

DO  J  -  Leadpt ( I ) + 1 .  Leadpt ( 1+ 1 )- 1 
Skip2  -  Skip2  +  1 
IF  (  MOD(Skip2,  2)  . NE .  0  )  THEN 
Count  -  Count  +  1 
X(  Level,  J  )  -  Mx(  Count  ) 

Y(  Level,  J  )  -  My(  Count  ) 

ENDIF 

ENDDO 

ENDIF 

ENDDO 


RETURN 


Listing  b.  Code  for  plotting  skeletal  fractal  mountains  (cont'd) 


This  subroutine  draws  the  latticeC Iter ) . 


SUBROUTINE  Draw(  X,  Y,  Iter  ) 

REAL *8  X( 0 : 6 ,  24768).  Y(0:6,  24768)  !  Array  of  endpoints 

INTEGER *4  Iter 

INTEGER *4  I.  J,  Depth.  Int.  Jump 
INTEGER*4  Leadpt(129) 

CHARACTER* 1  Lntyp  !  Draws  either  points(P)  or  vectors(V) 

100  FORMAT ( 1X.F10.3, 1X.F10.3, IX, A) 

Depth  -  2 
DO  I  =  1.  Iter 

Depth  =  Depth  +  2**(I-1) 

ENDDO 

Leadpt(l)  -  1 
Int  »  0 

DO  I  =  2,  ( Depth- 1) 

Int  -  Int  +  1 

Leadpt(I)  =  Leadpt(I-l)  +  Int 
ENDDO 

--  Drawing  line  segments  "slanting  left" 

DO  I  -  1.  ( Depth- 1) 

Lntyp  *=  ’  P ' 

WRITE( 10. 100)  x(  Iter,  ( Leadpt( 1+ 1 ) - 1 )  ). 

*  Y(  Iter.  (Leadpt(I+l)-l)  ),  Lntyp 
Jump  ~  0 

DO  J  -  1,  (Depth-I ) 

Jump  »(J+I-1)+  Jump 
Lntyp  =  'V' 

WRITE( 10. 100)  X(  Iter.  (  (Leadpt(I+l)-l)  +  Jump  )  ), 

*  Y(  Iter,  (  (Leadpt ( 1  +  1 ) - 1 )  +  Jump  )  ), 

*  Lntyp 
ENDDO 

ENDDO 

--  Drawing  line  segments  "slanting  right” 

DO  I  =  1.  ( Depth- 1) 

Lntyp  -  ' P ' 

WRITE( 10. 100)  X(  Iter.  Leadpt ( I )  ). 

*  Y(  Iter,  Leadpt(I)  ).  Lntyp 
Jump  -  0 

DO  J  *  1,  (Depth-I) 

Jump  -  (J  +  I)  +  Jump 
Lntyp  -  'V' 

WRITEC 10. 100)  X(  Iter.  (Leadpt(I)  +  Jump)  ). 

*  Y(  Iter,  (Leadpt(I)  +  Jump)  ),  Lntyp 
ENDDO 

ENDDO 

--  Drawing  "horizontal"  line  segments 
DO  I  -  2,  Depth 
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Listing  6.  Code  for  plotting  skeletal  fractal  mountains  (cont'd) 


WRITE ( 10 , 100)  X(  Iter,  Leadpt(I)  ), 

*  Y(  Iter,  Leadpt(I)  ),  Lntyp 

DO  J  -  (Leadpt(I)+l) ,  (Leadpt(I+l)-l) 

Lntyp  -  'V* 

WRITE(10, 100)  X(  Iter,  J  ), 

‘  Y(  Iter,  J  ),  Lntyp 

ENDDO 
ENDDO 

RETURN 

END 


5.  REMARKS  ON  APPLICATIONS 


Because  fractals  apparently  mimic  nature  so  well,  they  have  been  applied 
to  the  study  of  numerous  areas.  Chemists,  biologists,  physicists,  statiti- 
cians,  etc,  have  been  using  fractals  lately  to  model  behavior  in  their  parti¬ 
cular  fields.  Fractals  could  even  be  applied  to  digital  signal  processing,  as 
shown  by  the  following  argument. 

Fractals  are  quasi-self-repeating  images.  Therefore,  by  definition,  the 
image  seen  on  one  level  is  nearly  mimicked  by  an  enlargement  on  the  next. 
Moreover,  the  only  limitations  to  this  magnification  are  the  built-in  limita¬ 
tions  of  the  image-producing  machines  (the  fractal  itself  will  allow  any 
magnification) . 

Therefore,  aligned  in  the  proper  digital  fashion,  a  fractal  could  be  used 
as  a  communications  code.  The  digitized  fractal  image  would  act  like  a 
carrier  or  an  envelope  along  which  information  would  travel.  It  should  make  a 
good  code  for  the  following  two  reasons: 

Uni queness — A  fractal  image  has  its  own  unique  imprint.  Small  perturba¬ 
tions  in  the  fractal  can  produce  large  variation.  (This  can  be  demon¬ 
strated  by  varying  A  in  the  iteration  of  z2  +  A.)  Thus,  properly  chosen, 
a  fractal  can  produce  a  unique  digital  sequence. 

Stability  (and  instability) — Fractals  should  be  extremely  stable  with 
respect  to  noise  because  of  their  quasi-self-repeating  nature.  To  decide 
whether  or  not  a  bit  is  a  good  piece  of  information,  the  bit  could  be 
enlarged  at  a  lower  level.  There,  the  criterion  of  the  information's 
correctness  can  be  predetermined  by  how  far  it  is  from  the  fractal's  basic 
quasi-pattern.  Further  enlargements  only  improve  upon  the  fractal's 
estimate. 
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